On the effect of surfactant 
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Abstract: Substantial experimental, theoretical, as well as numerical effort has been 
invested to understand the effect of boundary slippage in microfluidic devices. However, 
even though such devices are becoming increasingly important in scientific, medical, 
and industrial applications, a satisfactory understanding of the phenomenon is still 
lacking. This is due to the extremely precise experiments needed to study the problem 
and the large number of tunable parameters in such systems. 

In this paper we apply a recently introduced algorithm to implement hydrophobic 
fluid-wall interactions in the lattice Boltzmann method. We find a possible explanation 
for some experiments observing a slip length depending on the flow velocity which is 
contradictory to many theoretical results and simulations. Our explanation is that a 
velocity dependent slip can be detected if the flow profile is not fully developed within 
the channel, but in a transient state. 

Further, we show a decrease of the measured slip length with increasing viscosity and 
demonstrate the effect of adding surfactant to a fluid flow in a hydrophobic microchan- 
nel. The addition of surfactant can shield the repulsive potential of hydrophobic walls, 
thus lowering the amount of slip with increasing surfactant concentration. 

Keywords: lattice Boltzmann, microflows, apparent slip 
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Introduction 



^yjicroflow devices are used for chemical, biological, or med- 
ial analysis techniques. Putting the "lab on a chip" al- 
J^ws to minimize the time needed for the analysis with only 
5«iall amounts of fluid. Also, such microdevces are more 

^mobile and allow a parallel treatment of multiple fluids. 

>*£>ther microflow systems are used as sensors and actors 
drjr devices like chemical reactors, cars, airplanes and inkjet 
printers. 

In these miniature apparatuses, a number of effects ap- 
pear which cannot easily be explained with our conven- 
tional physical understanding. A common example is the 
violation of the no-slip boundary condition. The no-slip 
boundary condition is one of the fundamental assumptions 
common in classical fluid mechanics, stating that the ve- 
locity of a fluid at a wall is equal to the velocity of the 
wall. For macroscopic applications no-slip is undoubted 
but during recent years a number of experiments found 
a violation of the no-slip boundary condition on a length 
scale of nanometers upto micrometers U 0). Numerous 
experiments (0; 0; H; i; E 0; 0; [1) utilize a modified 
atomic force microscope (AFM) with an oscillating col- 



loidal sphere at the tip of its cantilever to measure the 
force needed to displace the fluid between the colloidal 
sphere and a wall. From the detected force, the amount of 
wall-slippage can be estimated as described in Other 
authors like Tretheway and Meinhart apply particle image 
velocimetry (PIV) to observe the flow near the fluid-wall 
boundary directly to quantify wall slippage (flf3;[l~lh. How- 
ever, it is still an open question if the detected slip is a 
fundamental property or appears due to surface variations, 
uncertainties in the experimental setups, or the complex 
interactions between fluid and wall. 

Instead of the no-slip boundary condition, Navier in- 
troduced in 1823 a slip boundary condition where the 
transversal velocity near the wall v z [x = 0) is proportional 
to the shear rate and the so called slip length /3 $T* 



v z (x = 0)=(3^\ x=0 . 



(1) 



Here, the boundary is at x = 0. z is the flow direction 
and v z is the fluid velocity in flow direction, parallel to the 
wall. The slip length j3 can be interpreted as the distance 
between the wall and the virtual point inside the wall at 
which the extrapolated flow velocity would be zero. 



Due to the large number of tunable experimental param- 
eters like temperature, viscosity, flow velocity, pressure, or 
surface properties, as well as their individual dependencies 
on each other, it is not possible to cover all occuring phe- 
nomena in a single experiment. In fact, a change in the 
surface properties usually implies a different experimen- 
tal setup and a change of viscosity without varying the 
temperature is only possible by a replacement of the fluid. 
However, such strong interventions might also have an in- 
fluence on other parameters of the system. In computer 
simulations it is possible to vary a single parameter of the 
fluid, e.g., the viscosity or the density, without changing 
other parameters. This is important to improve our un- 
derstanding of the effects occuring in microffuidic systems 
and to further promote the design of such devices. 

In addition, computer simulations are able to study the 
properties of multiphase flows in microchannels with the 
individual fluid parameters and fluid-fluid interactions be- 
ing well defined. In particular, the influence of surfactant 
is of interest here. Surfactant molecules are often called 
amphiphiles and are comprised of a hydrophilic (water- 
loving) head group and a hydrophobic (oil-loving) tail. In 
a non-wetting microchannel filled with water, the surfac- 
tant molecules arrange at the interface between water and 
surface, thus shielding the hydrophobic repulsion of the 
wall. On the other hand, in a wetting channel an arrange- 
ment of surfactant molecules at the boundary causes the 
otherwise wetting wall to become hydrophobic. As a result 
an apparent slip occurs. 



particles" behaves as the corresponding real microscopic 
system. Due to this coarse-graining, the numerical effort 
is much smaller than for molecular dynamics simulations 
because the collision and propagation rules of the used 
"quasi particles" are much simpler than the ones of real 
particles. Therfore, much larger particle counts can be 
simulated for substantially longer times. An example for a 
mesoscopic simulation method is "stochastic rotation dy- 
namics" (SRD), which is sometimes called "multi particle 
collision dynamics" (MPCD). In a propagation step, every 
representative fluid particle is moved according to its veloc- 
ity to its new position. In the collision step, the simulation 
volume is split into boxes. In each box the velocity vectors 
of every single particle are rotated around the mean veloc- 
ity in a random manner, so that energy and momentum are 
conserved in every box (|2 ll l22i ). The method is efficient 



and is used when Brownian motion is required. Its disad- 
vantage is that thermal fluctuations cannot be switched off. 
"Dissipative particle dynamics" (DPD) also utilizes quasi 
particles which represent a set of molecules. The propaga- 
tion of such a collective quasi particle is implemented as 
in molecular dynamics but collisions are dissipative. This 
method is easy to implement in an existing MD simulation 
code but the computational costs are still very high. 

In this paper we use the lattice Boltzmann method, 
where one discretizes the Boltzmann kinetic equation 



d I 
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ry(x, v, t) = n 



(2) 



2 Simulation method 



The simulation method used to study microfluidic devices 
has to be choosen carefully. While Navier- Stokes solvers 
are able to cover most problems in fluid dynamics, they 
lack the possibility to include the influence of molecular 
interactions as needed to model boundary slip. Molecular 
dynamics simulations (MD) are the best choice to simulate 
the fluid- wall interaction, but the computer power today is 
not sufficient to simulate length and time scales necessary 
to achieve orders of magnitude which are relevant for ex- 
periments. However, boundary slip with a slip length f3 of 
the order of many molecular diameters a has been studied 
with molecular dynamics simulations by various authors 
(U; EI EE H El H; El)- They find increasing slip with 
decreasing liquid density and liquid-solid interactions as 
well as a decrease of slip with increasing pressure. How- 
ever, the maximum number of particles that can be simu- 
lated on today's most powerful supercomputers is about 20 
billion |20h. This corresponds to a volume of less then one 



/am , but the typical length scale of a microfluidic device 
is about 100/im. 

A mesoscopic model is able to govern a volume large 
enough to describe the flow properties and still holds in- 
formation about the molecular behavior. The term "meso- 
scopic" means that the trajectories of single molecules are 
not simulated in detail but a whole ensemble of "quasi 



on a lattice, rj indicates the probability to find a single 
particle with mass to, velocity v at the time t at position 
x. The derivatives represent simple propagation of a single 
particle in real and velocity space whereas the collision op- 
erator f2 takes into account molecular collisions in which a 
particle changes its momentum due to a collision with an- 
other particle. External forces F can be employed to imple- 
ment the effect of gravity or external fields. To represent 
the correct physics, the collision operator should conserve 
mass, momentum, and energy, and should be Gallilei in- 
variant. By performing a Chapman Enskog procedure, it 
can be shown that such a collision operator $7 reproduces 
the Navier-Stokes equation (23). In the lattice Boltzmann 
method the time t, the position x, and the velocity v are 
discretized. 

During the last years a number of attempts to simulate 
slip within the lattice Boltzmann method have been devel- 
oped. The most simple idea is to use a partial bounce back 
boundary condition |23h. While full bounce back leads to 
no slip, full reflection leads to full slip. Partial slip implies 
that a particle is reflected by the wall with the probabil- 
ity q, while it bounces back with probability (1 — q). As 
a result, a finite boundary slip can be observed. Nie et 
al. (|24l ) use a Knudsen-number dependent relaxation time 
in the vicinity of the wall to generate slippage in an ideal 
gas lattice Boltzmann model. 

Our attempt to gene rate slip involves a repulsive po- 
tential at the wall (|25l). This leads to a depletion zone 



near the wall with a reduced density resulting in an ap- 
parent slip at hydrophobic (non wetting) walls. Bcnzi et 
al. Hil) introduced a similar approach but the repulsion 
there decays exponentially while the potential we are us- 
ing only takes into account next neighbor lattice sites as 
described below. Our method is based on Shan and Chen's 
multiphase lattice Boltzmann method, i.e., the interaction 
between the surface and the fluid is simulated similar to 
the interactions between two fluid phases. This allows us 
to recycle our well tested parallel 3D multiphase lattice 
Boltzmann code, as it is presented in (|27l ) with only minor 
modifications. It is very advantaguous of our model that 
its parameters can be linked to experimentally available 
properties, namely the contact angle (f28h. 

The simulation method and our implementation of 
boundary conditions are described as follows. A multi- 
phase lattice Boltzmann system can be represented by a 
set of equations ((29I ) 

f? ?(x + c i ,t + l)-» ? ?(x,t)=nf, z = 0,l,...,6, (3) 

where ry"(x, t) is the single-particle distribution function, 
indicating the amount of species a with velocity Cj, at 
site x on a D-dimensional lattice of coordination number b 
(D3Q19 in our implementation), at time-step t. This is a 
discretized version of equation ([2]) without external forces 
F for a number of species a. For the collision operator fif 
we choose the Bhatnagar-Gross-Krook (BGK) form (|30f ) 



n? 



(4) 



where r Q is the mean collision time for component a and 
determines the kinematic viscosity 



2r Q - I 
6 ' 



(5) 



of the fluid. The system relaxes to an equilibrium distri- 
bution 77" eq which can be derived imposing restrictions on 
the microscopic processes, such as explicit mass and mo- 
mentum conservation for each species ( 31c 132c 1331 ) . In our 
implementation we choose for the equilibrium distribution 
function 



1 



(<vu) 2 

2rf 



2 c 2 



(c,-u) 

6c6 
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23 



(6) 



which is a polynomial expansion of the Maxwell distribu- 
tion, a are the velocity vectors pointing to neighbouring 
lattice sites. c s = 1/V3 is the speed of sound for the 
D3Q19 lattice. The macroscopic values can be derived 
from the single-particle distribution function r]f(x,i), i.e., 
the density rj a (x) of the species a at lattice site x is the 
sum over the distribution functions r)°'(x) for all lattice 
velocities cl 

^(x,t) =5>f(x,t). (7) 

i 

u a (x, t) is the macroscopic velocity of the fluid, denned as 



(8) 



Interactions between different fluid species are introduced 
following Shan and Chen as a mean field body force be- 



tween nearest neighbors (|34t |3 



F Q (x,i) 



5> a (x',i)(x'-x) , (9) 



(1 



-r,°(x,t)/r,o) is the S o-called effec- 



where ^ Q (x, t) 

tive mass with m being a reference density that is set to 
1 in our case (|34|). g a& is a force coupling constant, whose 
magnitude controls the strength of the interaction between 
component a and a. The dynamical effect of the force is 
realized in the BGK collision operator (T3J) by adding an 
increment Su a — r a F a frf 1 to the velocity u in the equilib- 
rium distribution function ^ . For the potential of the wall 
we attach the imaginary fluid "density" 7y wa11 to the first 
lattice site inside the wall. The only difference between 
^waii anc j an y th er fluid packages on the lattice rj a is that 
the fluid corresponding to r/ wa11 is only taken into account 
for in the collision step, but not in the propagation step. 
Therefore, we can adopt and the coupling constant 

5a,waii in order to tune the fluid-wall interaction, (fo^aii is 
kept at 0.08 throughout this paper if not mentioned oth- 
erwise and all values are reported in lattice units. Addi- 
tionally, we apply second order correct mid-grid bounce 
back boundary conditions between the fluid and the sur- 
face (23). Extending our model to a multi-relaxation time 
scheme would result in a more correct treatment of the 
boundaries, but the difference in the observed slip lengths 
is expected to be neglectable since interaction induced by 
the repulsive force between fluid and wall causes a sub- 
stantially larger effect. 

From molecular dynamics simulations it is known that 
the fluid-wall interactions causing a slip phenomenon usu- 
ally take place within a few molecular layers of the liquid 
along the boundary surface (|13l ; 



14; ha hj; 17; 18 



Our coarse-grained fluid wall interaction acts on the length 
scale of one lattice constant and does not take the molec- 
ular details into account. Therefore, our implementation 
is only able to reproduce an averaged effect of the inter- 
action and we cannot fully resolve the correct flow profile 
very close to the wall and below the resolution of a single 
lattice spacing. However, in order to understand the in- 
fluence of the hydrophobicity on experimentally observed 
apparent slip, it is fully sufficient to investigate the flow 
behavior on more macroscopic scales as they are accessi- 
ble for experimental investigation. Our method could be 
improved by a direct mapping of data obtained from MD 
simulations to our coupling constant <?a,wall allowing a di- 
rect comparison of the influence of liquid- wall interactions 
on the detected slip. This is a currently ongoing project 
and our results will be published elsewhere. 

Amphiphiles are introduced within the model as de- 
scribed in |36T ) and (37). An amphiphile usually possesses 
two different fragments, one being hydrophobic and one be- 
ing hydrophilic. The orientation of any amphiphile present 
at a lattice site x is represented by an average dipole vec- 
tor d(x, t). Its direction is allowed to vary continuously 
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Figure 1: Slip length (3 versus bulk velocity v (on a logarithmic scale), for different fluid- wall interactions 7y wa11 after 
a) t = 50000 and b) t = 15000 time steps. For better visibility, both figures share the same legend. The slip length 
is independent of the flow velocity after 50 000 timesteps and only depends on 7y wa11 (Fig. a)). After 15000 timesteps, 
however, the slip length starts at a critical minimum velocity and appears to rise with increasing v (Fig. b)). Even 
though the mean flow velocity has reached its final value already and the parabolic velocity profile is well developed, 
the system is still in a transient state at t = 15000 (see Fig. [3])) resulting in an eventually misleading measurement of 
p. All units are expressed in lattice units throughout this paper. 



and to keep the model as simple as possible no informa- 
tion is specified for velocities c.;. The surfactant density 
at a given site is given by an additional fluid species with 
density 77 sur ,that behaves as every other species a. The 
direction d(x, t) propagates with the fluid field according 
to 



3 Simulation setup 



if 



r (x,t+l)d(x,i+l) = ^r7f r (x-c^)d'(x-c^) (10) 



and during the collision step the direction d evolves to the 
equilibrium direction d eq similar to the BGK operator 

x ,/ x d(x,t) - d e «(x,i) 
d'x,f =dx,t - ^ ' V ' 1 11 

(d' indicates the direction after the collision step). The 
equilibrium distribution d eq ~ 4f-h is proportional to the 
so called color field or order parameter h which represents 
the distribution of the other species. It is defined as the 
weighted sum of the densities of all species 



(12) 



In our case (a = 2) we set the weights to e a — ±1, i.e., 
h corresponds to the density difference between the two 
species. 

The model has been used successfully to study spinodal 
decomposition ( 38l ; l39l) . binary and ternary amphiphilic 
fluids under shear (|40f ). the formation of mesophases (ill ; 
H; Hi; 0; 0; H3), and flow in porous media (H3). Of par- 
ticular relevance for the present paper is our first article on 
simulations of apparent slip in hydrophobic microchannels 
21 



The simulations in this work use a setup of two infinite 
planes separated by the distance 2d. We call the direction 
between the two planes x and if not stated otherwise 2d 
is set to 64 lattice sites. In y direction we apply periodic 
boundary conditions. Here, 8 lattice sites are sufficient to 
avoid finite size effects, since there is no propagation in 
this direction, z is the direction of the flow with our chan- 
nels being 512 lattice sites long. At the beginning of the 
simulation (t = 0) the fluid is at rest. We then apply a 
pressure gradient Vp in the z- direction to generate a pla- 
nar Poiseuille flow. Assuming Navier's boundary condition 
|T|) , the slip length f3 is measured by fitting the theoretical 
velocity profile, 



v z (x) 



J_dP 

2/x dz 



2dp] 



(13) 



in flow direction (v z ) at position x, to the simulated data 
via the slip length f3. We validate this approach by com- 
paring the measured mass flow rate J rjv(x)dx to the the- 
oretical mass flow without boundary slip and find a very 
good agreement. The pressure gradient is realized by 
a fixed inflow pressure (P(z = 0) = c 2 s ri(z = 0) = 0.3 
if not stated otherwise). At the outflow (z = z max ) 
we linearly extrapolate the density gradient by setting 
?7(z max ) = 277(z max - 1) - T](z max - 2) in order to simu- 
late infinite plates. Therefore, the body force regulates 
the velocity. The dynamic viscosity fj, as well as the pres- 
sure gradient needed to fit equation (fl3|) are obtained 
from our simulation data. 

In a previous work fldi) , we have shown that this model 
creates a larger slip (3 with stronger interaction, namely 



larger g a ,waii and larger jy wa11 . The relaxation time T a was 
kept constant at 1.0 in this study and the maximum avail- 
able slip length measured was 5.0 in lattice units. For 
stronger repulsive potentials, the density gradient at the 
fluid- wall interface becomes too large, causing the simula- 
tion to become unstable. At lower interactions the method 
is very stable and the slip length (3 is independent of the 
distance d between the two plates and therefore indepen- 
dent of the resolution. We have also shown that the slip de- 
creases with increasing pressure since the relative strength 
of the repulsive potential compared to the bulk pressure 
is weaker at high pressure. Therefore, the pressure reduc- 
tion near the wall is less in the high pressure case than in 
the low pressure one. Furthermore, we have demonstrated 
that (3 can be fitted with a semi analytical model based on 
a two viscosity model. 



4 Results 

We have studied the dependence of the slip length j3 on 
the flow velocity for a wide range of velocities of more than 
three decades as it can be seen in Fig. [T]a) and in (f25h - In 
the figure, we show data for different fluid-wall interactions 
< r] waU < 2.0 and flow velocities from 10~ 4 < v < 
10 _1 . Within this region we confirm the findings of many 
steady state experiments (|48t liih. i.e., that the slip length 
is independent of the flow velocity and only depends on the 
wettability of the channel walls. Experimentalists often 
present measurements for different shear rates S, which 
for Poiseuille flow are given by 



a_ 9u t 
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(14) 



Some dynamic experiments, however, find a shear rate 
dependent slip 



[50j;[5l|) 



These experiments often uti- 
lize a modified atomic force microscope as described in the 
introduction to detect boundary slippage. Since the slip 
length is found to be constant in our simulations after suf- 
ficiently long simulation times, we investigate the behavior 
of the slip during the transient, i.e., for simulation times 
i<f c with t c = L z /v being the self convection time. The 
flow that is initially at rest has not converged to its fi- 
nal steady state. The time development of the slip length 
could explain an apparent shear dependence as shown in 
Fig. [1] b), where j3 is plotted over the flow velocity for 
different fluid-wall interactions at t = 15000. Here, the de- 
tected (3 depends very strongly on the flow velocity. The 
figure shows a qualitative similarity to the data presented 
in (fHoh . namely there seems to be a critical shear rate at 
which the slip starts to increase very fast. However, this 
only holds during the transient as shown in Fig. [T]- in the 
steady state the slip is independent of the velocity. 

Fig. [2] depicts the time dependence of the measured slip 
length at constant ^ wal1 = 1.0 and for final flow velocities 
v = 0.7- 10~ 3 , 1.3 -10- 3 , and 4.0- 10- 3 . Since for t < 10000 
the expected parabolic velocity profile is not developed, we 
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Figure 2: Measured slip length (3 versus time t for different 
bulk velocities at constant 77 wa11 = 1.0. The difference be- 
tween the converged slip length and the slip length during 
the transient is greater for slower velocities. After the con- 
vection time t c = L z /v the slip is converged, but already 
for t > 50000 only small deviations from the final value 
can be observed. 



only plot our data for 10000 < t < 50000. It can be ob- 
served that the slip length develops to its final value for all 
three bulk velocities. However, the number of timesteps 
needed to achieve the steady state of (3 is dependent on v. 
The slip has reached its steady state after the convection 
time t c = L z /v, which is the time it takes for an individ- 
ual fluid element to be transported through the whole sys- 
tem. The slip converges with different rates depending on 
the flow velocity, but after 50000 timesteps the difference 
between the actual slip length and the converged one is 
neglectible already. This explains the fluctuations for very 
low velocities in Fig. . The determination of the correct 
slip length therefore can only be expected after sufficiently 
long simulation times. As can be seen from Fig. [3j it is not 
sufficient to just check if the velocity profile seems to have 
reached its final shape. Here, velocity profiles after 15000 
and 50000 timesteps are shown for a representative simula- 
tion run and 7y wa11 = 2.0. Even though the parabolic veloc- 
ity profile is already fully developed after 15000 timesteps, 
the measured slip length is (3 = 0.55 ± 7 ■ 10~ 3 only, while 
after 50000 timesteps /3 = 1.088±7 - 10~ 4 is obtained. The 
solid lines in Fig. correspond to a fit of the data with 
equation (fl~3)) . 

The kinematic viscosity v is another important param- 
eter in fluid dynamics. However, in experiments it can 
only be varied by changing the fluid itself and therefore 
it is inevitable to change other parameters too. Within 
the lattice Boltzmann method with BGK collision oper- 
ator ([5]), the kinematic viscosity of the fluid is given by 
([3]) and depends on the relaxation time r Q . Within the 
Shan-Chen model, a change of r Q also has an influence on 
the effect of the body force that enters the BGK opera- 
tor to model the fluid-fluid interactions. One might argue 
that this is not realistic since a change of viscosity does 
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Figure 3: The velocity profile v(x) for ry wa11 = 2.0 after 
t = 15000 and t = 50000 time steps. The lines are the 
parabolic fit with equation (fT3|) with a slip length of j3 = 
0.55 ± 7 • 10~ 3 at t = 15000. After 50000 time steps the 
slip length is significantly larger at f3 = 1.088 ± 7 • 10~ 4 . 
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Figure 4: Corrected slip length P(rj wa11 ) - /3(^ wa11 = 0.0) 
versus kinematic viscosity v a for 77 wa " = 0.5, 1.0, and 2.0. 
The slip length converges to as shown by the solid lines 
representing an exponential least squares fit of the data. 



not necessarily modify the fluid-fluid interactions between 
different species. Additionally, it is known that mid grid 
bounce back boundary conditions are second order correct 
while using the BGK collision operator, as it is used in 
this paper (|23l : |52|). For relaxation times r Q « 1 the error 
introduced due to the boundary condition is neglectible. 
However, we are interested in studying the dependence of 
boundary slippage on the fluid's viscosity. Therefore, we 
performed simulations with ^ wal1 = 0, i.e., without any 
fluid-wall repulsion, to estimate the effect of the error in- 
duced by the boundaries. For = 0, (3 should be zero 
as well, but we find the error of the slip length being pro- 
portional to (t") 2 - This behavior is expected by the theory 
of He et al. (|52r ) and can only be avoided by using a multi 
relaxation time approach. For 1 < r Q < 3 the numerical 



error is less than 5% of the slip length while for larger re- 
laxation times the error increases strongly so that the slip 
seems to increase. In order to reduce the influence of the 
error introduced by the single relaxation time method and 
the particular boundary conditions used, we subtract the 
slip length determined for ^ wal1 = from the measured 
(3 at ?7 wan > 0. The results are plotted in Fig. [H where 
we demonstrate a decreasing slip length with increasing 
viscosity for ?y wa11 = 0.5, 1.0, and 2.0. The data shown 
in Fig. [4] can be fitted exponentially as depicted by the 
solid lines and all three curves converge to zero for high 
viscosities. 
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Figure 5: Slip length /3 versus the concentration of surfac- 
tant in % for ?7 wa11 = 1.0. (3 is steadily decreasing with 
increasing the surfactant concentration from 0.64 down to 
0.19. The dashed line is given by a fit of the data with an 
exponential function. 

Since surfactant molecules consist of a hydrophobic and 
a hydrophilic part, they like to assemble at the interface be- 
tween a fluid and wetting or no n- wetting walls. As found 
by experimentalists, in a wetting microchannel, this can 
cause no slip to switch to partial slip |49t IBH ). In a non- 
wetting environment, the surfactant molecules can shield 
the hydrophobic repulsion of the surface H). We apply the 
amphiphilic lattice Boltzmann model as described earlier 
in this paper to model a fluid within a hydrophobic mi- 
crochannel that contains a surfactant concentration of up 
to 33%. The interaction parameters are choosen according 
to earlier works © EH; El EI H El; El), in such a way 
that they are not too strong to cause structuring effects 
in the flow, but strong enough to have an effect at the 
fluid-solid boundary. The total density inside our system 
r\ a + rj sur is kept fixed at 0.3. As initial condition the sys- 
tem is filled with a binary mixture of surfactant and fluid. 
The orientation d of the dipoles is choosen randomly. In 
Fig. 03 we plot the measured slip length for fluid- wall in- 
teractions determined by 77™ 11 = 0.5, 1.0 and 2.0 versus 
the concentration of surfactant. The symbols in Fig. [5] are 
given by the simulation data while the lines correspond to a 
fit with an exponential function. We find a strong decrease 
of the slip length with a higher surfactant concentration. 
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Figure 6: A typical profile of the surfactant concentra- 
tion in x direction, i.e., between the channel walls. Near 
the surface, the surfactant concentration is substantially 
higher (44%) than in the bulk (32%) since it is energeti- 
cally more favorable for the surfactant molecules to arrange 
at the flui-surface interface, thus shielding the repulsive po- 
tential of the hydrophobic channel wall. 

For all three values of ?7 wa11 , the measured slip lengths con- 
verge to the same value at high surfactant concentrations 
showing that at high concentrations the amount of surfac- 
tant that can assemble at the interface is saturating. 

In Fig.[6]we present a representative density profile of the 
surfactant for ?y wa11 = 1.0. The initial amphiphile concen- 
traton is set to 33% here. It can be seen that the concentra- 
tion at the first lattice site next to the surface increases to 
44%, while the bulk concentration stays constant at 32% 
- a value slightly lower than the initial 33%. This high 
concentration regime close to the boundary causes the hy- 
drophobic potential of the wall to be shielded and results 
in a decreasing slip. Our findings are consistent with ex- 
perimental results (f49l [|| [HT|) . 

Large amphiphilic molecules or polymer brushes show a 
shear dependent slip (f53h since they have to align with the 
shear forces acting on them. The higher the shear force, 
the more they are rotated causing the effect of shielding 
the hydrophobicity to be reduced. Since in our model the 
amphiphiles are point-like, we cannot expect to observe 
any shear rate dependence of (3. 



5 Conclusion 

In conclusion we have presented three-dimensional multi- 
phase lattice Boltzmann simulations which govern a wide 
range of slip phenomena. After demonstrating the validity 
of our model, we presented studies of the dependence of 
the boundary slip on the flow velocity. While the slip is 
independent of the velocity if the system is in the steady 
state, we find an apparent velocity dependence during early 
times of the simulation. For small numbers of timesteps, 
the parabolic velocity profile is already well developed, but 



due to the system being in a transient state, the detected 
slip is not correct. This is an important finding for exper- 
imental setups since to the best of our knowledge only dy- 
namic experiments find a velocity dependence, while static 
experiments confirm the slip lengths being independent of 
the flow velocity. Our findings are in good agreement with 
most non dynamic experiments |l|; 0) and MD simulations 

(USES)- 

For experimentalists it is a major effort to change the 
viscosity of the fluid without changing any other parame- 
ters of their setup. In computer simulations, however, this 
can be done easily. In our simulations we found a decrease 
of the detected slip with increasing viscosity. 

With a simple dipole model we were able to simulate 
the shielding of the repulsive potential between hydropho- 
bic walls and a fluid if surfactant is present in the solu- 
tion, i.e., the slip length decreases with increasing surfac- 
tant concentration. However, we were not able to show a 
shear dependence as it is seen in experiments with polymer 
chains. In a future work, we plan to extend our simula- 
tions to govern larger molecules which can be affected by 
a shear flow. Then, we hope to be able to study the shear 
rate dependence of boundary slippage. 
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